function [mu,V,df] = mlstudent(Y,df)

% This matlab function returns maximum-likelihood estimates of mu, sigma, and degrees 
% of freedom using data following a multivariate t-distribution with unknown df.

%%% Inputs:
% Y: a TxN matrix of data
% df: initial estimate of degrees of freedom (default = 10)

%%%O utputs:
% mu: ML estimate of mu
% V: ML estimate of covariance matrix V
% df: ML estimate of degrees of freedom

maxiter = 300;
OPTIONS = optimset('Display','off');
[T,p] = size(Y);
if nargin<2
   df = 10;	%Initial estimate
end
Z = ones(T,1)*(df-2)/df;
w = Z/sum(Z);	
for iter=1:maxiter
    mu = w'*Y;
    Y1 = Y-ones(T,1)*mu;
    psia = (Y1'*(Z*ones(1,p).*Y1))/T;
    delta = sum(Y1.*(Y1*inv(psia)),2);
    df = fzero(@f2,df,OPTIONS);
    Z = (df+p)./(df+delta);	%Conditional expectation of Z
    w = Z/sum(Z);
    if all(abs(w'*Y-mu)<1e-10)
       break
    end
end
V = psia*(df/(df-2));
if iter==maxiter
   fprintf(' Warning!  Maximum number of iterations reached.\n')
   fprintf(' Degrees of freedom = %8.1f\n',df)
end

function y = f2(x)
   if x<=0
      y = 10000-x;%Make sure x stays away from negative
   else
      Z = (x+p)./(x+delta);
      y = psi(0.5*(x+p))-psi(0.5*x)+log(x/(x+p))+mean(log(Z)-Z)+1;
   end
end

end
